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ABSTRACT 

Kilometer-sized moonlets in Saturn's A ring create S-shaped wakes called "propellers" 
in surrounding material. The Cassini spacecraft has tracked the motions of propellers 
for several years and finds that they deviate from Kcplerian orbits with constant 
semimajor axes. The inferred orbital migration is known to switch sign. We show 
using a statistical test that the time series of orbital longitudes of the propeller Bleriot 
is consistent with that of a time-integrated Gaussian random walk. That is, Bleriot's 
observed migration pattern is consistent with being stochastic. We further show, using 
a combination of analytic estimates and collisional iV-body simulations, that stochastic 
migration of the right magnitude to explain the Cassini observations can be driven 
by encounters with ring particles 10-20 m in radius. That the local ring mass is 
concentrated in decameter-sized particles is supported on independent grounds by 
occultation analyses. 

Key words: diffusion - methods: statistical - methods: numerical - planets and 
satellites: rings - planet-disc interactions - celestial mechanics 



1 INTRODUCTION 

Propellers are disturbances in Sa turn's rings caused b y 
moonlets 0.1-1 km in radius (|Tiscareno et al. I l2006t ). 
The moonlets gravitationally repel material to either 
side of their orbits, creating partial gaps that diffuse 



shut via inter-particle c ollisions (ISpahn fc Sremcevidl2000l; 
Sremcevic et alj |2002| ; ISeifi et alj 120051 : iLewis fc Stewart] 



2009h . 



Even with the Cassini spacecraft's resolving power, 
the moonlets themselves are too small to be detected di- 
rectly. Their existence and sizes are inferred from the 
larger S-shaped wakes resembling propelle rs that they leave 
behind, on scales of several Hill radii |Seifi et al. 120051 : 
iTiscareno et alj|2010j . hereafter T10). 

Intriguingly, multi-epoch Cassini observations of a num- 
ber of propellers reveal that propellers deviate from strictly 
Keplerian orbits: the orbital longitudes A of a propeller drift 
away from the values expected for an orbit of fixed semi- 
major axis (T10). Longitude residuals AA range from 0.01- 
0.31 degrees over At — 1.3-4.3 years; see Table 1 of T10. By 
far the most extensive data exist for the propeller dubbed 
Bleriot, whose longitudes have been measured 89 times at 
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sporadic intervals over 4.2 years. Bleriot's longitude resid- 
uals versus time are shown in Figure [TJ reproduced from 
T10. The measurements of non-Keplerian motion represent 
the first direct evidence that moons embedded in rings ex- 
hibit orbital evolution. Propellers thus provide a test-bed 
for studying satellite-disk in teractions, in particular orbi tal 
migration (see, for example. iGoldreich fc Tremainell 19821 ). 

To establish the order of magnitude of the implied mi- 
gration, we convert longitude residual AA to the change Aa 
in the moonlet's semimajor axis: 
AA 



Aa 



flAt 



30 



AA 
0.1 deg 



2yr 
At 



1.3 x 10 5 km 



E/2 



(1) 



where Q is the mean motion (orbital frequency). This in- 
ferred radial deviation should not be confused with the az- 
imuthal deviation, which is measured directly from observa- 
tions to be on the order of aAA ~ 300 km. The semimajor 
axis change Aa observed to date is smaller than the moon- 
let's expected physical size. 

Two classes of theories have emerged to explain 
the non-Keplerian motion. In one, propeller-moonlets li- 
brate within a resona nce established by co-orbital material 
jPan fc Chiang||2010l ). In the other, propeller-moonlets are 
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Figure 1. Longitude residuals (deviations from a fixed circular orbit) of the propeller Bleriot as seen in 89 Cassini images obtained 
over 4.2 years and as reported by T10. For display purposes, we subtracted a linear trend to place the first and last data points at zero 
longitude residual. 



torqued stocha stically by density fluctua t ions in surroundin g 
ring material (|Rein fc Papaloizoul l20ld ; ICrida et al.l |20T0| ) . 
Both types of theories in their current forms are not with- 
out problems. Libration amplitudes within the proposed co- 
orbital resonance damp to zero when account is made of 
the two-way feedbac k between the moonlet and the ring 
|Pan fc Chiandfeoilj ). So far, the stochastic migration hy- 
pothesis has focused on den sity fluctua t ions driven by self- 
gravit y ("self-gravity wakes";[§alo 1995). lRein fc Papaloizoul 
|201fJ ) simulated such wakes, and found that they could 
cause small moonlets, about 25-50 m in radius, to ran- 
dom walk by distances comparable to those cited above. 
Unfortunately, most propellers, including Bleriot, are an or- 
der of magnitude larger in size (T10), and thus are not ex- 
pected to be accelerated significantly by self-gravity wakes. 
This shortcoming of self-gravity wakes can be remedied 
by in creasing the ring surface density (e.g., ICrida et al.l 
I2010D, but evidence for such large surface densities is lacking 
|Colwell et afl|200gh . 

It may seem surprising that even the most basic issue of 
whether the non-Keplerian motion is deterministic or ran- 
dom is controversial. A glance at the longitude time series 
in Figure [1] suggests that Bleriot 's motion is smooth or even 
sinusoidal, with a period of ~3.7 yr. Nevertheless, the time 
series could in fact reflect pure noise. The confusion arises 
because orbital longitude is a time-integrated quantity: 

AA = J AQ dt (2) 

where AQ is the difference between the moonlet's instan- 
taneous mean motion and that of a fixed reference orbit. 
The time integration smooths over fluctuations in AQ; AX 
is obviously differentiable. The integration also introduces 
correlations between data points even when AQ itself rep- 



resents uncorrelated noise: AA(t) depends on the full time 
history of perturbations up to time t. Both effects conspire 
to hide any underlying stochasticity. 

In this paper, we apply a test, well-known among statis- 
ticians but less so among astronomers, that identifies inte- 
grated random walks for the special case where the under- 
lying random walk is Gaussian. In effect, the test "undoes" 
the correlations introduced by the integration to determine 
whether the integrand AQ is consistent with a Gaussian 
random walk. This test, which we call the "diagonalization 
test" for reasons explained in [J2] is especially useful because 
it can be applied to data that — as is the case for real-life 
propellers — are unevenly sampled. Q 

The plan of this paper is as follows. In [JU we describe 
the diagonalization test and apply it to Bleriot. We find that 
Bleriot passes the test — that its behavior is consistent with 
that of an integrated Gaussian random walk. In fJ3]we de- 
scribe how such a random walk can be driven by Poisson 
fluctuations in the encounter rate between large ring parti- 
cles and a given propeller moonlet (an effect distinct from 
self-gravity wakes). A summary is given in Sj4] 



1 If the sampling were uniform, then we could in principle just 
histogram second differences of the sequence of points to check 
consistency with an integrated Gaussian random walk. Taking 
differences between adjacent points would produce the Gaussian 
random walk observed on a uniform grid of times. Taking differ- 
ences again would produce a sequence of independent, identically 
distributed Gaussian random variables, with each one being the 
sum of some common number of the Gaussian steps in the under- 
lying random walk, taken over a set of times disjoint from those 
leading to the other sums. 
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2 DIAGONALIZING THE NOISE 



timesteps, n < m, is the expected value of their product: 



We derive and explain the rationale behind the diagonaliza- 
tion test in £12.11 check the test against some case examples 
in £12.21 and apply the test to Bleriot in £12.31 



2.1 Decorrelating the Integrated Gaussian 
Random Walk 

We wish to check if a given time series — for example, 
the Bleriot longitude residuals — is consistent with a time- 
integrated random walk whose individual steps are inde- 
pendent identically-distributed (IID) Gaussian random vari- 
ables. We can think of this time-integrated random walk 
as correlated Gaussian noise in which all the correlations 
arise from the time integration. If all such correlations were 
to be eliminated from our given time series, we could then 
compare what remains to a family of IID Gaussian random 
variables and thereby test the Gaussian random walk hy- 
pothesis. Here we describe one such decorrelation method. 
The mathematical ideas behind the m ethod are well known 
(see, for example. iMardia et al.|[l979l ). Indeed, the method 
overlaps c onceptually with recent treatments of pulsar tim- 
ing noise (|Coles et al.ll201ll) . 

We first calculate the covariance between any two ele- 
ments of a time series produced by a time-integrated Gaus- 
sian random walk. Although our presentation focuses on the 
specific case of a time series of orbital longitudes, the un- 
derlying algorithm is general. We use A for quantities that 
evolve with time and 8 for quantities associated with an in- 
dividual step in the random walk. 

Consider a body whose semimajor axis a undergoes a 
random walk, taking one step in the walk per timestep. After 
n ^ 1 timesteps of constant length 8t, the total change in 
semimajor axis is 



n-l 

Aa{nSt) = ]T & 

j=0 



(3) 



where the are Gaussian random variables for which 
= an d E^j^k] = (6a) 2 Sjk- Here E denotes the ex- 
pected value, and Sjk is the Kronecker delta. 

The longitude residual AA is the time integration of 
Aadn/da: 



on 

AX(n8t) = -^2 —Aa(iSt)5t 



2a 



3QSt \ \ -* 
~^2a~^^ 3 

i=l j=0 



3mt 

2a 



(4) 
(6) 



3=0 



where, as before, St is the time interval between steps. 

The covariance between values of AA after n and m 



E[AX(nSt)AX(mSt)} 



3=0 



my* 



( snstSa 

1 2a 



^2^3<3(n-j)(m-j) 

3=0 
2 n-l 

3=0 



(ZQMSa^ 
\ 2a 



^(1 — n 2 + 3m + 3nm) 



(7) 

(8) 
(0) 

(10) 

(11) 



Line (£9|) follows because AX(nSt), AX(mSt) are 
two snapshots of the same integrated random walk 
{AX(St), AX(2St), AX(n5t), AX(mSt), ...}. Between 
timesteps 1 and n, the histories of AX(n8t) and AX(m5t) 
coincide exactly — they are 100% correlated — so for 
summation indices j, k < n, in effect i = j. After timestep 
n, the Gaussian variables . . . ,£ m -i contributing to the 
further history of AX(mSt) are completely independent of 
those that contributed to AX(nSt), so the determinants 
of the motion after timestep n do not contribute to 
E[AX(nSt)AX(mSt)]. 

Given a list of observations at times {tk : 1 ^ k ^ 
K}, and choosing St to be the characteristic time between 
changes in semimajor axis, we can use the above with 
n — tk./St and m = ti/St to calculate the entries of the cor- 
responding (positive definite, symmetric) covariance matrix 
(Efc^i^fsgK- If we write AA* for the column vector with fc th 
entry AX(t k ), then E is just the K x K matrix E[aX A~X t ]. 
By the spectral theorem, there is an orthogonal matrix U 
and a diagonal matrix Z with positive diagonal entries such 
that E = U ZU T . To eliminate the correlations in A A due to 
the time integration, we simply use this covariance matrix 
to define a suitable linear transformation of the time-series 
vector AA: 

diagonalized (decorrelated) residuals f 

= UZ- 1/2 U T Ai, (12) 

where Z^ 1 ^ 2 Z~ 1//2 = Z -1 , the inverse of Z. The covariance 
matrix of such a column vector f is 

E[rr T ] = UZ~ 1/2 U T E[Ai Ak T ]UZ~ 1/2 U T (13) 

= UZ~ 1/2 U T T,UZ~ 1/2 U T (14) 

= UZ~ 1/2 U T UZU T UZ- 1/2 U T (15) 

= /, (16) 

the K xK identity matrix. Because linear transformations of 
Gaussian rando m vectors are also G aussian (see, for exam- 
ple, Chapter 3 of lMardia et al.iri979T ). the entries in r are IID 
Gaussian random variables with common expected value 
and common standard deviation 1 — this fact i s effect ively 
the content of Corollary 3.2.1.1 of lMardia et all (|l979h . 

Note that the covariance matrix E is in practice un- 
known, as only the AA are observed. However, E is of the 
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form E = c 2 E, where c = and E fc£ = |(l-n 2 + 3m + 

3nm) with n = tk/St and m = til St. The matrix E is known, 
but the value of c is unknown (it depends, for example, on 
the unknown step size Sa) and so c must be estimated from 
data. We do this as follows. 

We have E = UZU T , where U and Z = c~ 2 Z are read- 
ily computed. Set 

f = UZ' 1/2 U T A\ = cr. (17) 

By the argument above, f is a vector of independent Gaus- 
sians with common mean and common standard devia- 
tion c. We estimate c using the standard estimator c = 

\J ~k 5^fc=i ^fc- T ne vec t° r c _1 f should then be close to r 
when K is not too small and hence should have entries that 
are approximately Gaussian with common expected value 
and common standard deviation 1. We can check this by 
plotting the empirical cumulative distribution function of 
the entries of c" 1 ? against the cumulative distribution func- 
tion of a Gaussian distribution with expected value and 
standard deviation 1. 

In essence, we have re-expressed the vector of measure- 
ments of an integrated Gaussian random walk at the times 
{tk} in a new basis so that the coefficients with respect to 
the new basis are independent and identically distributed 
— hence our term "diagonalization test" . We expect that if 
the {A\(tt)} arise from a two- fold time integration of indi- 
vidual IID Gaussian kicks, then the entries of c _1 f ~ f will 
also be distributed as IID Gaussians. 

Thus, the diagonalized residuals c _1 f provide a conve- 
nient negative test of whether a time series is consistent with 
an integrated Gaussian random walk. If the entries of the 
vector c _1 f derived from a given observation vector A A do 
not follow a Gaussian distribution reasonably closely, then 
the time series A^ cannot result directly from an integrated 
Gaussian random walk. Conversely, if the entries of c _1 f are 
approximately Gaussian, then the observations {AA(t^)} are 
consistent with (but do not uniquely demand) an integrated 
Gaussian random walk. 

The procedure outlined above is predicated on the as- 
sumption that the longitude residuals {AX(tk)} are observed 
without measurement uncertainty, so that the covariance 
matrix of the observations is some multiple of E. We will 
examine the effects of measurement uncertainty in £|2 . 2 1 and 
H2.3I The diagonalization procedure also requires that the 
timestep St be less than or equal to the time interval between 
any two of the observation times {tk}, so that the columns 
of E are independent. We show in H2.3l that any sufficiently 
small St that satisfies this condition will result in essentially 
the same covariance matrix E (and hence essentially the 
same probability model for A^), provided that Sa is mod- 
ified so that the diffusivity D — (Sa) 2 /St is kept constant. 
Moreover, we show that any sufficiently small choice of St 
will give essentially the same value of c _1 f and essentially 
the same diagonalization test. However, for Bleriot there is 
a natural, physically motivated choice of timestep, St = n _1 
(see 

Finally, for easy reference we give a short summary of 
the diagonalization test algorithm: 

(i) The input data required are a time series of values 



{AA(ffc)} and the corresponding list of times {tk}; the goal 
is to check if the {AX(tk)} are consistent with an integrated 
Gaussian random walk. 

(ii) Assuming a constant St smaller than the shortest in- 
terval ife+i — tk, use Efc£ = ^(1 — n 2 + 3m + 3nm) with 
n — tk/St and m = ti/St to compute the matrix E. 

(iii) Compute U and Z from E = UZU T using standard 
matrix decomposition procedures. 

(iv) Compute f from U, Z, and the {AA} using Eq. 1171 

(v) Estimate c using c = <J -L XIa-Li t\, and check that 
the distribution of c _1 f has a Gaussian shape and a standard 
deviation of about 1. 

A simple Mathematica implementation of the diagonaliza- 
tion test is available in an online supplement to this paper. 



2.2 Gaussian Walks vs. Levy Flights 

As a simple check of the diagonalization test, we apply it to 
a simulated integrated Gaussian random walk0 The param- 
eters of our simulation are motivated by Bleriot. We take the 
time interval between steps to be St = Q^ 1 ~ 8842 s and 
integrate the walk for 4.3 years. Semimajor axis changes, 
or "kicks" Sa to the moonlet, are generated as Gaussian 
random variables with standard deviation 1 m. We time- 
integrate 5a(t) once to get the cumulative semimajor axis 
evolution Aa(t), and twice to get the associated longitude 
variations AA(t). The simulated longitudes are then sampled 
at the times of the Bleriot observations. When the Bleriot 
data are appropriately binned (see £|2 .3 P , there are 41 obser- 
vation times. 

Figure [2] shows the results of the diagonalization test 
when applied to our simulated time series containing 41 
points. The diagonalized residuals c~ 1 {fk} appear close to 
Gaussian. In our calculation of c we dropped 4 extreme val- 
ues among the {?k} because including these outliers skewed 
c so as to be clearly inconsistent with the vast majority of the 
{ffc}. Aside from these few outliers, which actually fall out- 
side the range of the right-hand panel in Figure [2] the agree- 
ment with a Gaussian is satisfactory. We find c ~ 1.2c where 
c is computed using the input parameters of the simulated 
random walk, which we also consider satisfactory agreement. 

To better calibrate what we mean by "satisfactory," we 
also apply the diagonalization test to data that does not de- 
rive from an integrated Gaussian random walk. We generate 
instead a time-integrated Levy flight, where the steps in the 
underlying random walk are drawn from a power law. Specif- 
ically, we assume that the probability of getting a kick of size 
Sa or larger scales as |5a|~ 2,/5 . This power law distribution 
describes kicks excited by perturbers that are sparsely dis- 
tributed over an annulus much wider than the moon let 's Hill 
sphere radius (see, for example. ICollins fc SarilfeoOrj . and ref- 
erences therein). Just as in the Gaussian experiment above, 
we time-integrate the kicks twice to get the associated lon- 
gitude variations, and we sample the longitudes at the 41 



2 The simulations described here in [[2] should not be confused 
with the 3D shearing box simulations of [|3] 
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Figure 2. Checking the diagonalization test with a simulated integrated Gaussian random walk. A kick drawn from a Gaussian distribu- 
tion of width 1 m in semimajor axis is applied every St = = 8842 s. The left panel shows the simulated longitude residuals sampled 
at the Bleriot observation times (binned). For display purposes, we subtracted a linear trend to place the first and last data points at 
zero longitude residual. The right panel shows the results of the diagonalization test. The sorted diagonalized residuals (filled circles) are 
a reasonable match for Gaussian random variables (solid line). The 4 outliers dropped in calculating c all lie outside the range shown in 
this plot. 




time (yrs) diogonolized residual 



Figure 3. Checking the diagonalization test with a simulated integrated Levy flight: same as Figure [2] except that here the kick 
magnitudes \8a\ have a power-law cumulative distribution oc | <5a. | 2 / 5 _ sor t e( j diagonalized residuals (filled circles in right panel) do 
not conform to a Gaussian distribution (solid line). In this case, 3 outliers were excluded in the calculation of c, and all fall outside the 
range shown in the right-hand panel. 



binned Bleriot observation times. We then perform the di- 
agonalization test on the samples. As Figure [3] shows, the 
diagonalized residuals are distinctly non-Gaussian in shape. 

In a world free of measurement uncertainty, we could 
use c to constrain the product of 5a /a and QSt and thus 
obtain a joint constraint on the moonlet's diffusivity D = 
(8a) 2 /St. Unfortunately, we have found by direct experiment 
that c is quite sensitive to measurement uncertainty. For ex- 
ample, randomly altering the longitudes shown in the left 
panel of Figure [2] by ~2%, an amount comparable to the 
uncertainties reported for the observed Bleriot longitudes, 
gives c an order of magnitude larger than that computed 
using the unaltered longitudes. A similar result is obtained 



for the Levy flight experiment in Figure [3] Fortunately, the 
shapes of the diagonalized residual distributions still yield 
the same qualitative answers: with measurement uncertain- 
ties included, the simulated integrated Gaussian random 
walk still passes the diagonalization test, and the simulated 
integrated Levy flight still fails the test. Thus we remain 
confident that, for the parameters of the problem at hand, 
as long as relative measurement uncertainties remain at the 
level of a few percent — as they seem to be for the actual 
data of Bleriot — the diagonalized residuals can still be used 
to give a "yes-or-no" answer to the question of whether the 
input data are consistent with an integrated Gaussian ran- 
dom walk. 
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2.3 Bleriot 

We apply the diagonalization test to Bleriot. In mapping a 
given time tk to an integer n = tk/St, we take the width St of 
each time bin to equal the dynamical time QT 1 = 8842 s, for 
the physical reason that each encounter between the moonlet 
and a ring particle which changes the moonlet's semimajor 
axis takes a dynamical time to complete (see i)3.1[) . We could 
take time bins of larger width, but that would reduce the 
number of points in our input time series. The precise choice 
of St is, in any case, not crucial because the random walk 
is approximately scale invariant: the distribution of A^ de- 
pends primarily on the combination of St and Sa through 
the diffusion coefficient D = (8a) 2 /St. To see this, note from 
Eq. [11] that the covariance matrix E of A A satisfies 



Efci = 



3QSt8a\ t k 
2a ) 65t 



st +6 st +6 stst 



for tk ^ ti. Hence, the probability model for AA is, to first 
order, unaffected by a change in St, provided that St is suffi- 
ciently small and Sa is also changed so as to keep the diffu- 
sivity D constant. The same calculation shows that if St is 
changed to hSt for some constant h, then E is, to first order, 
changed to /i~ 3 E. Then f and c become respectively /i 3//2 f 
and h 3 ' 2 c, so that c _1 f is, to first order, unchanged. Conse- 
quently, the diagonalization test is effectively invariant with 
respect to a change in the choice of St. 

Many points in Bleriot 's published time series are sepa- 
rated by less than St = Q~ , as multiple images were taken 
in short succession. The longitudes that fall into the same 
time bin are averaged into one number. Measurement un- 
certainties are added in quadrature. Binned this way, there 
are 41 "independent" longitude measurements, as shown in 
Figure H (left panel) 

The matrix operations in Eq. [12] are applied to the 
binned longitude residuals {AA(^)} to compute the diago- 
nalized residuals c~ {ffc}. In this case, 3 outlier values were 
dropped in the computation of c. Figure [4] shows the dis- 
tribution of diagonalized residuals. It is close to Gaussian; 
compare with Figures [2] and [3] 

We have explored the sensitivity of these results to mea- 
surement uncertainties in Bleriot's data. Ten different real- 
izations of Bleriot's longitude time series were generated by 
randomly selecting points within the error bars shown in the 
left panel of Figure [4] In all cases, the diagonalized residuals 
resembled those shown in the right panel of Figure 2] 

Recall our experiments in H2. 21 with simulated random 



3 Within some bins, the original longitudes differ by more than 
their published error bars (see, e.g., the cluster of five points near 
year 2005.38 in Figure 4b of T10). The published uncertainties are 
probably underestimated, especially in those cases where longi- 
tudes are referenced to ring features such as the Encke gap edge 
instead of to stars. We have verified with M. Tiscareno (2012, 
personal communication) that our binning procedure is justified 
given the quality of the data. 



walks, where we found that the magnitude of the diagonal- 
ized residuals was sensitive to measurement uncertainty; in 
particular, the (unnormalized) diagonalized residuals were 
much larger than expected when measurement errors were 
added to the random walks. Indeed the same effect seems 
to manifest here with the actual data for Bleriot. The c 
value for Bleriot is an order of magnitude larger than the 
c value for the simulated integrated Gaussian walk even 
though their longitude residuals are of comparable magni- 
tudes. Thus we have no reliable estimator of the true value 
of c, and thus no constraint on the moonlet's diffusivity D 
from the results of the diagonalization test alone. 

Although we cannot measure the moonlet's diffusivity 
from the diagonalization test, we still have the original lon- 
gitude time series of Bleriot, in addition to a smattering 
of longitude data for other propellers (see Table 1 of T10). 
Taken at face value, these data indicate that moonlet dif- 
fusivities must be such as to generate "typical" longitude 
deviations of ~0. 1-0.3 deg over timescales of ~l-2 yr. In 
the next section, we explain how such random walks can 
physically arise. 



3 STOCHASTIC MIGRATION DUE TO LOCAL 
FLUCTUATIONS IN SURFACE DENSITY 

A moonlet can undergo a Gaussian random walk because 
of Gaussian fluctuations in the number of particles encoun- 
tered at Hill sphere separations. The fluctuations of interest 
to us do not depend on ring self-gravity. Fluctuations still 
occur because of frequent collisions between ring particles 
that randomize their positions in the time intervals between 
encounters with the moonlet. We sketch this process ana- 
lytically ( H3.ip . test and calibrate our order-of-magnitudc 
scaling relations with numerical iV-body simulations f i]3.2[) . 
and apply our theory of stochastic migration to Bleriot and 
other propellers ( i|3.3| l. Many of the ideas in this section have 

Chiandl2006l ; 



been treated previously ( e.g., Murray-Clay 
iRein fc Papaloizoull2010l : ICrida et alJl20ld ) 



|), but we present 



them here afresh for clarity and convenience. 



3.1 Analytic Description of Gaussian Stochastic 
Migration 

Consider a moonlet of radius i? m00 n at semimajor axis a, 
embedded in a ring composed of particles each of radius r 
and mass m. The moonlet and ring particles are assumed 
individually to have a bulk density pi,. The surface mass 
density of the ring is E, and the local orbital frequency is 
Q. Ring particles shear by the moonlet and gravitationally 
perturb it. Particles inside the moonlet's orbit tend to kick 
the moonlet onto a larger orbit, while particles outside the 
moonlet's orbit tend to push the moonlet inward. 

Random fluctuations in the rate of particles encoun- 
tered cause the moonlet's semimajor axis to change stochas- 
tically. At radial separations on the order of x between a 
collection of ring particles and the moonlet, the relative Ke- 
plerian shearing velocity is ~ Qx. The duration of an en- 
counter is 8t onc ~ x/(Qx) ~ independent of x. The 
number of particles passing conjunction (to either side of 
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Figure 4. Results of the diagonalization test for Bleriot: same as Figurc[2]but for the actual data for Bleriot rather than for simulation 
results. The data are binned into 41 points such that the time interval between bins is at least f2 _1 = 8842 s long. The diagonalized 
residuals are reasonably close to Gaussian-distributed; compare with Figures [2] and [3] As in previous figures, the 3 outliers excluded from 
the calculation of c fall outside the range shown in the right panel. 



the moonlet's orbit) per St enc should follow a Poisson distri- 
bution with mean N enc ~ T,x 2 /m and width V N cnc . 

The randomly varying excess number of particles — ly- 
ing to either side with equal probability — is responsible 
for net changes in the moonlet's semimajor axis. We call the 
mass of this excess group of particles, tallied every encounter 
time 5t e nc, the "fluctuation mass" mfl uct . If JV cn c S> 1, we 
can treat mfl uct as approximately Gaussian-distributed with 
mean zero and width ^m^/ N anc : positive/negative signs cor- 
respond to excess groups of particles passing outside/inside 
the moonlet's orbit. An encounter with a fluctuation mass 
occurs once every Stenc, and each encounter changes the 
moonlet's velocity by 8v ~ (Gmj uct /i 2 ) x 8t cnc , with G 
the gravitational constant. The largest fluctuations arise 
fro m particles within severa l Hill s phere radii of the moon- 
let (|Murrav-Clav fe Chianel l200rj ; ICrida et all |201Ch : x ~ 



an rms distance 



Rmu = a (^ L pi,7?moon/"iSaturn) 1/3 - For such encounters, the 
fractional change in the moonlet's semimajor axis is of order 
the fractional change in its velocity: 8a/ a ~ 8v/(Qa), with 
equal probability of either signQ 

Putting all of the above together, we find that every 
8t cnc ~ SI -1 time, the moonlet randomly walks in semimajor 
axis by a step of root-mean-square size 



8a • 



WT-fluct 
WSaturn \ ^?Hill 

' 300 m 



0.15 



Rn 



) (lOm) 



40 g/cm 2 



3/2 

!/ 2 /l / 3\ 1/6 

' 1 g/cm x 



f>b 



(18) 



Over the course of the Cassini observations analyzed by T10, 
a propeller-moonlet will random walk in semimajor axis by 



4 This relation holds for x ~ Rum, but not for x 3> ^?mn- 



Aa(t) ~ SaVnt ~ 10 



300 m 

R moon 



(lOm) 



3/2 



t 



1/2 



(19) 

where we have used S = 40 g cm 2 and pt — 1 g cm 3 . 
Comparison with Eq. [T] shows that this is of the right order 
of magnitude to explain the observed non-Keplerian motions 
of propellers. In the next section we employ A^-body simu- 
lations to test our scaling relations and measure more accu- 
rately the moonlet's diffusivity. That is, we will revise our 
crudely estimated coefficients in Eqs. [18] and [19] in accord 
with the numerical simulations. 



3.2 Numerical Simulations of Gaussian Stochastic 
Migration 

We perform shearing box simulations of a moonlet randomly 
perturbed by nearby ring particles. We use the freely avail - 
able collisional A'-body code REBOUND l|Rein fc Liul I2012T I 
with shear periodic boun dary conditions and the symplectic 
epicycle integrator (SEI. lRein fc Tremainell2011r ). 

The simulatio ns are similar to th e test-particle simula- 
tions performed bv lCrida et ail (|2010h but include particle- 
particle a nd particle-moonlet collis ions. Unlike the simu- 
lation s of iRein fc Papaloizoul (|2010l ') and iLewis fc Stewart! 
|2009[ ). ours do not explicitly include self-gravity. As ex- 
plained in fJTJ self-gravity wakes probably make only a small 
contribution to the observed longitude residuals of propellers 
as large as Bleriot. The mean self-gravitational field does 
enhance the vertical freque ncy f2 z of epicyclic motion (e.g., 
IWisdom fc Tremaind 119881 ); we mock up this effect in the 
code by setting Q z — 3.6S1. All simulations are performed 
with Q = 1.131 • 10~ 4 s 
axis of a = 130000 km. 



1 corresponding to a semimajor 
The timestep was chosen to be 



dt = 10" 3 2tt/Q. 



Simulation parameters are listed in Table [T] The ring 
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Figure 5. Distribution of kicks (changes in scmimajor axis) felt 
by our simulated 100-m moonlet embedded in a ring with surface 
density £ = 40 g cm -2 , maximum particle size r max = 10 m, min- 
imum particle size r m ; n = 2.5 m, and power-law index q = 3 for 
the differential size distribution of particles. Kicks are computed 
every time interval Q~ ; their cumulative distribution (solid line) 
conforms closely to a Gaussian (open circles). The typical kick 
size, 5a ~ 1.5 m, is consistent with our order-of-magnitude esti- 
mate in Eq. \TS\ 



particles and moonlet are assumed to have a bulk density 
pb = 0.4 g/cm 3 . Ring particles are assumed to follow a dif- 
ferential size distribution dN/dr oc r~ q from r m i n to r max . 
The slope q is fixed at 3, which places most of the mass 
in the largest ring particles. Some of our chosen particle 
size parameters are compatible with occultation and imag- 
ing observations (Cuzzi et al. 2009; T10). Others were chosen 
only to provide a large enough dynamic range to probe how 
stochasticity scales with the particle radius r. 

Ring particles are initialized with zero random velocity 
(i.e., their initial velocity is determined purely by Keplerian 
shear). Once a ring particle exits an azimuthal boundary 
of the simulation domain, it re-enters the domain on the 
opposite side at a randomized radial location (semimajor 
axis), with zero random velocity. Thus the number of ring 
particles N in the box remains constant. The dimensions of 
the box are L x , L y , and L z in the radial, azimuthal, and 
vertical directions, respectively. In nearly all cases, L x x Ly 
covers ~27 x 135 moonlet Hill radii, while L z is chosen large 
enough so that no particle ever reaches a vertical boundary. 

Figure [5] shows the distribution of semimajor axis 
changes or "kicks" 5a, evaluated every time interval St — 
l/fi, for simulation 9. The empirical distribution of kicks 
is close to Gaussian, confirming our physical description of 
stochastic migration in £13. II Furthermore, the characteris- 
tic value of 5a (defined as the la half- width of the Gaussian 
distribution) is 1.5 m, which agrees to order-of-magnitude 
with the prediction of Eq. [18] for a 100-m moonlet. In fact, 
the simulated characteristic value for 5a is about 3 times 
larger than predicted by our back-of-the-envelope estimate. 
The enhanced fluctuations revealed by the numerical simula- 
tions strengthen the case for perturbations from the largest 
(decameter-sized) ring particles ( £|3.3[1 as the main cause of 
the propellers' observed non-Keplerian motions. We spec- 



ulate that the factor of 3 might arise from two effects in 
our numerical simulations that are omitted from our sim- 
ple analysis in H3.ll ring particles on horseshoe orbits, and 
direct collisions between ring par ticles and the moonlet 
dRein fc Papaloizoul l201fj : see also iMurrav-Clav fc Chiang! 
120061 who show that encounters with horseshoe librators en- 
hance stochasticity by a factor of order unity). 

Figure [6] shows how the moonlet 's longitude residuals 
evolve in simulation 17, whose parameters are the same as 
those of simulation 9 but is run for 5 years. In order to avoid 
impractically long runtimes, simulation 17 uses a smaller box 
size than simulation 9. However, we expect the results to be 
insensitive to the change in box size: the most important 
interactions are with particles passing within ~i?mn of the 
moonlet ( H3.1[l and the smaller box is still several Rum in 
size. We applied the diagonalization test to simulation 17, 
sampling the longitude residuals at the same set of 41 times 
that characterize Bleriot's binned data. The diagonalized 
residuals, shown in Figure [6] verify that the behavior of the 
simulated moonlet is consistent with that of an integrated 
Gaussian random walk. 

Figure[7]shows the characteristic value of 5a versus both 
particle radius r and moonlet radius -Rmoon- We compute 
this characteristic value by plotting a cumulative Gaussian 
distribution against the kicks sorted from smallest to largest, 
performing a least-squares fit to a line, and varying the width 
of the Gaussian until the slope of the fitted line is 1. We take 
the Gaussian width so derived to be the characteristic 5a. 
The scalings expected from Eq. [18] are approximately con- 
firmed. The agreement for r is better than for -R m00 n, but we 
consider both acceptable. Using our numerical simulations 
to normalize our analytic scalings, we calibrate Eqs. 1181 and 
[T!J]into more accurate forms: 




(21) 

3.3 Implications for Bleriot 

Observations of Bleriot require Aa ~ 30 m over a period 
of ~2 years. We can use Eq. 1211 which is calibrated using 
numerical simulations, to estimate how big the surrounding 
ring particles must be to reproduce these observed parame- 
ters. 

The radius of Bleriot's moonlet is thought to lie in the 
range i? moon = 300-1200 m (see Figure 2 of T10). If we 
adopt i? m oon = 700 m, then a particle size of r ~ 18 m 
would satisfy the observations assuming a bulk density of 1 
g/cm 3 for all bodies. That is, the 1-cr excursion in semimajor 
axis for a 700-m moonlet is Aa ~ 30 m over At = 2 yr when 
r = 18 m. If the observed Aa ~ 30 m actually represents 
a 2-cr excursion, then the required particle size decreases to 
r = 11 m. 
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Table 1. Parameters of simulations using REBOUND, a shearing box code for colliding particles. The duration of each simulation is given 
in units of 1/Q. See text for description. 
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Figure 6. Results of the diagonalization test for simulation 17: same as Figure [2] but for simulation 17 sampled at the 41 binned 
observation times of Bleriot. Here 4 outliers were excluded from the calculation of c, and all lie outside the range shown in the right-hand 
panel. The longitude residuals shown in the left panel are simulation data sampled at times corresponding to the Bleriot observation 
times. The residuals' distribution is close to Gaussian. 



In fact, ground- and space-based occultation data of 
the outer A ring independently indicate that the bulk of 
the ring mass resides in particles of size r = 10-20 m (for 
a summary of what is known about pa rticle size distribu - 
tions based on occultation analysis, see ICuzzi et al.ll2009h . 
For the region just outside the Encke gap which contains 
Bleriot and the other giant propellers, fits to Voyager ob- 
servations yiel d a maximum par ticle size r max = 8.9 m 
and q = 3.03 (jZebker et al.lfl985h . An analysis of ground- 
based occultation data giv es r max = 20 m and q — 2.9 
IIFrench fc Nicholson! |200Ch . We conclude that stochastic 
gravitational interactions between propeller moonlets and 
the largest nearby particles in the outer A ring can readily 
reproduce longitude residuals like those observed for Bleriot. 



4 SUMMARY 

Whether the migration patterns of propellers arise from 
a deterministic or random process is not obvious just by 
looking at their longitude residuals. The difficulty arises be- 
cause longitude residuals are time-integrated quantities. The 
time integration smooths out semimajor axis variations that 
could be noisy, and introduces correlations between data at 
a given time and all earlier times. 

The "diagonalization test" removes correlations intro- 
duced by time integration of a Gaussian random process. 
It tests whether a given time series is compatible with an 
integrated Gaussian random walk. We have applied the diag- 
onalization test to the longitude time series of the propeller 
Bleriot and found that it passes the test. Bleriot's behavior 
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Figure 7. Characteristic 8a value, or "kick size", for simulations 1—11 (left panel) and 12—16 (right panel). Simulations 1-11 show how 
varying the maximum particle size affects the typical 8a. The best-fit power law shown has index 1.47. This agrees well with the predicted 
slope of 1.5 from Eq. 1181 The fit uses only data from simulations 6—11: in simulations 1—5, the largest ring particles are so rare that 
A'enc < 1. Simulations 12-16 show how varying the moonlet size affects the typical 8a. The best-fit power law shown has slope —0.74; 
given the scatter in the data, we consider this acceptable agreement with Eq. 1181 which predicts a slope of —1. 



is consistent with that of an integrated Gaussian random 
walk. 

By combining simple analytic scaling relations with nu- 
merical iV-body simulations, we also showed that moonlets 
as large as Bleriot, having radii of ~700 m, could exhibit 
longitude residuals on the order of 0.1 deg over 2 years, 
when embedded in a ring of surface density 40 g/cm 2 — 
provided the bulk of the mass of the ring is contained in 
particles 10-20 m in radius. Such ring properties are inferred 
on in dependent grounds by occultation analysis (|Cuzzi et al.l 
120091 ). The perturbations exerted by large ring particles on 
propeller-moonlets are stochastic, caused by Poisson fluc- 
tuations in the number of ring particles that shear by the 
moonlet on Hill sphere scales. 

The picture of stochastic m igration that we support i s 
similar to that first proposed bv lRein fc Papaloizoul (|2010t ). 
except that the primary contributors to stochasticity are 
decameter-sized particles, not self-gravity wakes. We have 
shown by direct iV-body simulation (e.g., Figure|6]) that par- 
ticle size distributions that place most of the ring mass in 
decameter sizes can reproduce longitude residuals like those 
observed. As the Cassini spacecraft emerges from the ring 
plane in 2012 and resumes observations of propellers, we look 
forward to measurements of longitude time series for other 
propellers in addition to Bleriot — and to more accurate pro- 
tocols for making longitude measurements by improvements 
to the matrix describing the orientation of the ISS camera 
with respect to the spacecraft. These new and more accu- 
rate data can also be subjected to the diagonalization test, 
and used to test our prediction that longitude residuals scale 
with moonlet size as A A oc -Rmoon- 
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